-- Finite slab of thickness L of magneto-optic material.
-- Assume the permittivity tensor in the receiving halfspace is
--   [  e   ie' 0 ]
--   [ -ie' e   0 ]
--   [  0   0   e ]
-- Then solving the wave equation gives two circular polarization modes:
--   E = (1, i)/sqrt(2), km = omega sqrt(e-ie)
--   E = (1,-i)/sqrt(2), kp = omega sqrt(e+ie)
-- We assume then that inside the slab there are waves of the form
--   Ex,y(z) = a(1,i)/sqrt(2) exp(i km z) + b(1,-i)/sqrt(2) exp(i kp z) + c(1,i)/sqrt(2) exp(-i km z) + d(1,-i)/sqrt(2) exp(-i kp z)
-- For normal incidence, we match boundary conditions:
--   k0 = omega
--   (E0+Erx, Ery) = a(1,i)/sqrt(2) + b(1,-i)/sqrt(2) + c(1,i)/sqrt(2) + d(1,-i)/sqrt(2)
--   k0(E0-Erx, -Ery) = a km(1,i)/sqrt(2) + b kp(1,-i)/sqrt(2) - c km(1,i)/sqrt(2) - d kp(1,-i)/sqrt(2)
--   a(1,i)/sqrt(2) e^{i km L} + b(1,-i)/sqrt(2) e^{i kp L} + c(1,i)/sqrt(2) e^{i km L} + d(1,-i)/sqrt(2) e^{i kp L} = (Etx, Ety)
--   a km(1,i)/sqrt(2) e^{i km d} + b kp(1,-i)/sqrt(2) e^{i kp d} - c km(1,i)/sqrt(2) e^{i km d} - d kp(1,-i)/sqrt(2) e^{i kp d} = k0(Etx, Ety)
-- Solving these gives
--[[
sol2 = Solve[{
     {E0 + Erx, Ery} == 
      a {1, I}/Sqrt[2] + b {1, -I}/Sqrt[2] + c {1, I}/Sqrt[2] + 
       d {1, -I}/Sqrt[2],
     {E0 - Erx, -Ery} == 
      a rm {1, I}/Sqrt[2] + b rp {1, -I}/Sqrt[2] - 
       c rm {1, I}/Sqrt[2] - d rp {1, -I}/Sqrt[2],
     a {1, I}/Sqrt[2] Exp[I km L] + b {1, -I}/Sqrt[2] Exp[I kp L] + 
       c {1, I}/Sqrt[2] Exp[-I km L] + 
       d {1, -I}/Sqrt[2] Exp[-I kp L] == {Etx, Ety},
     a rm {1, I}/Sqrt[2] Exp[I km L] + 
       b rp {1, -I}/Sqrt[2] Exp[I kp L] - 
       c rm {1, I}/Sqrt[2] Exp[-I km L] - 
       d rp {1, -I}/Sqrt[2] Exp[-I kp L] == {Etx, Ety}
     }, {Erx, Ery, a, b, Etx, Ety, c, d}] // Simplify;


{{Erx -> (E0 ((-1 + E^(2 I km L)) (-1 + E^(2 I kp L)) k0^4 - (-1 + E^(
           2 I km L)) (-1 + E^(2 I kp L)) km^2 kp^2 + 
        k0 km kp ((-1 + E^(2 I km L)) (1 + E^(2 I kp L)) km + (1 + E^(
              2 I km L)) (-1 + E^(2 I kp L)) kp) + 
        k0^3 (-(1 + E^(2 I km L)) (-1 + E^(2 I kp L)) km - (-1 + E^(
              2 I km L)) (1 + E^(2 I kp L)) kp)))/(((-1 + E^(
           2 I km L)) k0^2 - 
        2 (1 + E^(2 I km L)) k0 km + (-1 + E^(2 I km L)) km^2) ((-1 + 
           E^(2 I kp L)) k0^2 - 
        2 (1 + E^(2 I kp L)) k0 kp + (-1 + E^(2 I kp L)) kp^2)), 
  Ery -> (I E0 k0 (km kp ((-1 + E^(2 I km L)) (1 + E^(
              2 I kp L)) km - (1 + E^(2 I km L)) (-1 + E^(
              2 I kp L)) kp) + 
        k0^2 ((1 + E^(2 I km L)) (-1 + E^(2 I kp L)) km - (-1 + E^(
              2 I km L)) (1 + E^(2 I kp L)) kp) - (-1 + E^(
           2 I km L)) (-1 + E^(2 I kp L)) k0 (km^2 - kp^2)))/(((-1 + 
           E^(2 I km L)) k0^2 - 
        2 (1 + E^(2 I km L)) k0 km + (-1 + E^(2 I km L)) km^2) ((-1 + 
           E^(2 I kp L)) k0^2 - 
        2 (1 + E^(2 I kp L)) k0 kp + (-1 + E^(2 I kp L)) kp^2)), 
  Etx -> (
   2 E0 k0 (-E^(I (km + 2 kp) L) km (k0 - kp)^2 - 
      E^(I (2 km + kp) L) (k0 - km)^2 kp + E^(I kp L) (k0 + km)^2 kp +
       E^(I km L) km (k0 + kp)^2))/(((-1 + E^(2 I km L)) k0^2 - 
      2 (1 + E^(2 I km L)) k0 km + (-1 + E^(2 I km L)) km^2) ((-1 + 
         E^(2 I kp L)) k0^2 - 
      2 (1 + E^(2 I kp L)) k0 kp + (-1 + E^(2 I kp L)) kp^2)), 
  a -> -((Sqrt[2] E0 k0 (k0 + km))/((-1 + E^(2 I km L)) k0^2 - 
     2 (1 + E^(2 I km L)) k0 km + (-1 + E^(2 I km L)) km^2)), 
  b -> -((Sqrt[2] E0 k0 (k0 + kp))/((-1 + E^(2 I kp L)) k0^2 - 
     2 (1 + E^(2 I kp L)) k0 kp + (-1 + E^(2 I kp L)) kp^2)), 
  Ety -> (2 I E0 k0 (-E^(I (km + 2 kp) L) km (k0 - kp)^2 + 
      E^(I (2 km + kp) L) (k0 - km)^2 kp - E^(I kp L) (k0 + km)^2 kp +
       E^(I km L) km (k0 + kp)^2))/(((-1 + E^(2 I km L)) k0^2 - 
      2 (1 + E^(2 I km L)) k0 km + (-1 + E^(2 I km L)) km^2) ((-1 + 
         E^(2 I kp L)) k0^2 - 
      2 (1 + E^(2 I kp L)) k0 kp + (-1 + E^(2 I kp L)) kp^2)), 
  c -> (Sqrt[2] E^(2 I km L)
     E0 k0 (k0 - km))/((-1 + E^(2 I km L)) k0^2 - 
    2 (1 + E^(2 I km L)) k0 km + (-1 + E^(2 I km L)) km^2), 
  d -> (Sqrt[2] E^(2 I kp L)
     E0 k0 (k0 - kp))/((-1 + E^(2 I kp L)) k0^2 - 
    2 (1 + E^(2 I kp L)) k0 kp + (-1 + E^(2 I kp L)) kp^2)}}

Code generated by:


CForm[Simplify[ComplexExpand[Im[Simplify[First[c /. sol2 /. {E0 -> 1}]]]]]]

]]




S = S4.NewSimulation()
S:SetLattice({1,0}, {0,1})
S:SetNumG(1)
f = 0.4
S:SetFrequency(f)

eps  = 10
epsp = 3
L = 1
S:AddMaterial("Dielectric", {
	{eps, 0}, {0, epsp}, {0, 0},
	{0, -epsp}, {eps, 0}, {0, 0},
	{0, 0}, {0, 0}, {eps, 0}
	})
S:AddMaterial("Vacuum", {1,0})

S:AddLayer('AirAbove', 0, 'Vacuum')
S:AddLayer('Stuff', 1, 'Dielectric')
S:AddLayerCopy('AirBelow', 0, 'AirAbove')

S:SetExcitationPlanewave(
	{0,0}, -- incidence angles (spherical coordinates: phi in [0,180], theta in [0,360])
	{0,0}, -- s-polarization amplitude and phase (in degrees)
	{1,0}) -- p-polarization amplitude and phase

k0 = 2*math.pi*f
km = k0*math.sqrt(eps-epsp)
kp = k0*math.sqrt(eps+epsp)
rm = math.sqrt(eps-epsp)
rp = math.sqrt(eps+epsp)
cosm = math.cos(2*km*L)
cosp = math.cos(2*kp*L)
sinm = math.sin(2*km*L)
sinp = math.sin(2*kp*L)

function Cos(x)
	return math.cos(x)
end
function Sin(x)
	return math.sin(x)
end
function Power(a,b)
	return a^b
end
function Sqrt(x)
	return math.sqrt(x)
end

for z = -5,5,0.01 do
	Exr,Eyr = S:GetEField({0,0,z})
	
	Ex = 0
	Ey = 0
	if z < 0 then
		xr = (2 + 6*Power(rm,2) + 6*Power(rp,2) - 6*Power(rm,4)*Power(rp,2) - 6*Power(rm,2)*Power(rp,4) - 2*Power(rm,4)*Power(rp,4) +      2*(-1 + Power(rm,2))*(1 + 3*(1 + Power(rm,2))*Power(rp,2) + Power(rm,2)*Power(rp,4))*Cos(2*km*L) - (-1 + Power(rm,2))*(-1 + Power(rp,2))*(-1 + Power(rm,2)*Power(rp,2))*Cos(2*(km - kp)*L) - 2*Cos(2*kp*L) -      6*Power(rm,2)*Cos(2*kp*L) + 2*Power(rp,2)*Cos(2*kp*L) - 2*Power(rm,4)*Power(rp,2)*Cos(2*kp*L) + 6*Power(rm,2)*Power(rp,4)*Cos(2*kp*L) + 2*Power(rm,4)*Power(rp,4)*Cos(2*kp*L) + Cos(2*(km + kp)*L) -      Power(rm,2)*Cos(2*(km + kp)*L) - Power(rp,2)*Cos(2*(km + kp)*L) + Power(rm,4)*Power(rp,2)*Cos(2*(km + kp)*L) + Power(rm,2)*Power(rp,4)*Cos(2*(km + kp)*L) - Power(rm,4)*Power(rp,4)*Cos(2*(km + kp)*L))/(2.*(1 + 6*Power(rm,2) + Power(rm,4) - Power(-1 + Power(rm,2),2)*Cos(2*km*L))*(1 + 6*Power(rp,2) + Power(rp,4) - Power(-1 + Power(rp,2),2)*Cos(2*kp*L)))
		xi = (2*rm*(-1 + Power(rm,2))*(1 + 6*Power(rp,2) + Power(rp,4))*Sin(2*km*L) + (-1 + Power(rp,2))*      ((-1 + Power(rm,2))*(rm - rp + Power(rm,2)*rp - rm*Power(rp,2))*Sin(2*(km - kp)*L) + 2*(1 + 6*Power(rm,2) + Power(rm,4))*rp*Sin(2*kp*L) -         (-1 + Power(rm,2))*(-rp + Power(rm,2)*rp + rm*(-1 + Power(rp,2)))*Sin(2*(km + kp)*L)))/   (2.*(1 + 6*Power(rm,2) + Power(rm,4) - Power(-1 + Power(rm,2),2)*Cos(2*km*L))*(1 + 6*Power(rp,2) + Power(rp,4) - Power(-1 + Power(rp,2),2)*Cos(2*kp*L)))
		Ex = xr * math.cos(-k0*z) - xi*math.sin(-k0*z) + math.cos(k0*z)
		yr = (-2*rm*(-1 + Power(rm,2))*(1 + 6*Power(rp,2) + Power(rp,4))*Sin(2*km*L) + (-1 + Power(rp,2))*      ((-1 + Power(rm,2))*(-rp + Power(rm,2)*rp + rm*(-1 + Power(rp,2)))*Sin(2*(km - kp)*L) + 2*(1 + 6*Power(rm,2) + Power(rm,4))*rp*Sin(2*kp*L) -         (-1 + Power(rm,2))*(rm - rp + Power(rm,2)*rp - rm*Power(rp,2))*Sin(2*(km + kp)*L)))/   (2.*(1 + 6*Power(rm,2) + Power(rm,4) - Power(-1 + Power(rm,2),2)*Cos(2*km*L))*(1 + 6*Power(rp,2) + Power(rp,4) - Power(-1 + Power(rp,2),2)*Cos(2*kp*L)))
		yi = (-6*Power(rm,2) - 2*Power(rm,4) + 6*Power(rp,2) - 6*Power(rm,4)*Power(rp,2) + 2*Power(rp,4) + 6*Power(rm,2)*Power(rp,4) +      2*(-1 + Power(rm,2))*(Power(rm,2) + 3*Power(rp,2) + 3*Power(rm,2)*Power(rp,2) + Power(rp,4))*Cos(2*km*L) + (-1 + Power(rm,2))*(Power(rm,2) - Power(rp,2))*(-1 + Power(rp,2))*Cos(2*(km - kp)*L) +      6*Power(rm,2)*Cos(2*kp*L) + 2*Power(rm,4)*Cos(2*kp*L) + 2*Power(rp,2)*Cos(2*kp*L) - 2*Power(rm,4)*Power(rp,2)*Cos(2*kp*L) - 2*Power(rp,4)*Cos(2*kp*L) - 6*Power(rm,2)*Power(rp,4)*Cos(2*kp*L) +      Power(rm,2)*Cos(2*(km + kp)*L) - Power(rm,4)*Cos(2*(km + kp)*L) - Power(rp,2)*Cos(2*(km + kp)*L) + Power(rm,4)*Power(rp,2)*Cos(2*(km + kp)*L) + Power(rp,4)*Cos(2*(km + kp)*L) -      Power(rm,2)*Power(rp,4)*Cos(2*(km + kp)*L))/(2.*(1 + 6*Power(rm,2) + Power(rm,4) - Power(-1 + Power(rm,2),2)*Cos(2*km*L))*(1 + 6*Power(rp,2) + Power(rp,4) - Power(-1 + Power(rp,2),2)*Cos(2*kp*L)))
		Ey = yr * math.cos(-k0*z) - yi * math.sin(-k0*z)
	elseif z > L then
		zz = z
		z = z-L
		xr = (-2*(-2*Power(rm,2)*(1 + 6*Power(rp,2) + Power(rp,4))*Cos(km*L) + Power(rm,2)*Power(-1 + Power(rp,2),2)*Cos((km - 2*kp)*L) + Power(rp,2)*Cos((2*km - kp)*L) - 2*Power(rm,2)*Power(rp,2)*Cos((2*km - kp)*L) +        Power(rm,4)*Power(rp,2)*Cos((2*km - kp)*L) - 2*Power(rp,2)*Cos(kp*L) - 12*Power(rm,2)*Power(rp,2)*Cos(kp*L) - 2*Power(rm,4)*Power(rp,2)*Cos(kp*L) + Power(rp,2)*Cos((2*km + kp)*L) -        2*Power(rm,2)*Power(rp,2)*Cos((2*km + kp)*L) + Power(rm,4)*Power(rp,2)*Cos((2*km + kp)*L) + Power(rm,2)*Cos((km + 2*kp)*L) - 2*Power(rm,2)*Power(rp,2)*Cos((km + 2*kp)*L) +        Power(rm,2)*Power(rp,4)*Cos((km + 2*kp)*L)))/((1 + 6*Power(rm,2) + Power(rm,4) - Power(-1 + Power(rm,2),2)*Cos(2*km*L))*(1 + 6*Power(rp,2) + Power(rp,4) - Power(-1 + Power(rp,2),2)*Cos(2*kp*L)))
		xi = (2*rm*(1 + Power(rm,2))*(1 + 6*Power(rp,2) + Power(rp,4))*Sin(km*L) - rm*(1 + Power(rm,2))*Power(-1 + Power(rp,2),2)*Sin((km - 2*kp)*L) + rp*Sin((2*km - kp)*L) - 2*Power(rm,2)*rp*Sin((2*km - kp)*L) +      Power(rm,4)*rp*Sin((2*km - kp)*L) + Power(rp,3)*Sin((2*km - kp)*L) - 2*Power(rm,2)*Power(rp,3)*Sin((2*km - kp)*L) + Power(rm,4)*Power(rp,3)*Sin((2*km - kp)*L) + 2*rp*Sin(kp*L) +      12*Power(rm,2)*rp*Sin(kp*L) + 2*Power(rm,4)*rp*Sin(kp*L) + 2*Power(rp,3)*Sin(kp*L) + 12*Power(rm,2)*Power(rp,3)*Sin(kp*L) + 2*Power(rm,4)*Power(rp,3)*Sin(kp*L) - rp*Sin((2*km + kp)*L) +      2*Power(rm,2)*rp*Sin((2*km + kp)*L) - Power(rm,4)*rp*Sin((2*km + kp)*L) - Power(rp,3)*Sin((2*km + kp)*L) + 2*Power(rm,2)*Power(rp,3)*Sin((2*km + kp)*L) - Power(rm,4)*Power(rp,3)*Sin((2*km + kp)*L) -      rm*Sin((km + 2*kp)*L) - Power(rm,3)*Sin((km + 2*kp)*L) + 2*rm*Power(rp,2)*Sin((km + 2*kp)*L) + 2*Power(rm,3)*Power(rp,2)*Sin((km + 2*kp)*L) - rm*Power(rp,4)*Sin((km + 2*kp)*L) -      Power(rm,3)*Power(rp,4)*Sin((km + 2*kp)*L))/((1 + 6*Power(rm,2) + Power(rm,4) - Power(-1 + Power(rm,2),2)*Cos(2*km*L))*(1 + 6*Power(rp,2) + Power(rp,4) - Power(-1 + Power(rp,2),2)*Cos(2*kp*L)))
		Ex = xr * math.cos(k0*z) - xi * math.sin(k0*z)
		yr = (-2*rm*(1 + Power(rm,2))*(1 + 6*Power(rp,2) + Power(rp,4))*Sin(km*L) + rm*(1 + Power(rm,2))*Power(-1 + Power(rp,2),2)*Sin((km - 2*kp)*L) + rp*Sin((2*km - kp)*L) - 2*Power(rm,2)*rp*Sin((2*km - kp)*L) +      Power(rm,4)*rp*Sin((2*km - kp)*L) + Power(rp,3)*Sin((2*km - kp)*L) - 2*Power(rm,2)*Power(rp,3)*Sin((2*km - kp)*L) + Power(rm,4)*Power(rp,3)*Sin((2*km - kp)*L) + 2*rp*Sin(kp*L) +      12*Power(rm,2)*rp*Sin(kp*L) + 2*Power(rm,4)*rp*Sin(kp*L) + 2*Power(rp,3)*Sin(kp*L) + 12*Power(rm,2)*Power(rp,3)*Sin(kp*L) + 2*Power(rm,4)*Power(rp,3)*Sin(kp*L) - rp*Sin((2*km + kp)*L) +      2*Power(rm,2)*rp*Sin((2*km + kp)*L) - Power(rm,4)*rp*Sin((2*km + kp)*L) - Power(rp,3)*Sin((2*km + kp)*L) + 2*Power(rm,2)*Power(rp,3)*Sin((2*km + kp)*L) - Power(rm,4)*Power(rp,3)*Sin((2*km + kp)*L) +      rm*Sin((km + 2*kp)*L) + Power(rm,3)*Sin((km + 2*kp)*L) - 2*rm*Power(rp,2)*Sin((km + 2*kp)*L) - 2*Power(rm,3)*Power(rp,2)*Sin((km + 2*kp)*L) + rm*Power(rp,4)*Sin((km + 2*kp)*L) +      Power(rm,3)*Power(rp,4)*Sin((km + 2*kp)*L))/((1 + 6*Power(rm,2) + Power(rm,4) - Power(-1 + Power(rm,2),2)*Cos(2*km*L))*(1 + 6*Power(rp,2) + Power(rp,4) - Power(-1 + Power(rp,2),2)*Cos(2*kp*L)))
		yi = (2*(2*Power(rm,2)*(1 + 6*Power(rp,2) + Power(rp,4))*Cos(km*L) - Power(rm,2)*Power(-1 + Power(rp,2),2)*Cos((km - 2*kp)*L) + Power(rp,2)*Cos((2*km - kp)*L) - 2*Power(rm,2)*Power(rp,2)*Cos((2*km - kp)*L) +        Power(rm,4)*Power(rp,2)*Cos((2*km - kp)*L) - 2*Power(rp,2)*Cos(kp*L) - 12*Power(rm,2)*Power(rp,2)*Cos(kp*L) - 2*Power(rm,4)*Power(rp,2)*Cos(kp*L) + Power(rp,2)*Cos((2*km + kp)*L) -        2*Power(rm,2)*Power(rp,2)*Cos((2*km + kp)*L) + Power(rm,4)*Power(rp,2)*Cos((2*km + kp)*L) - Power(rm,2)*Cos((km + 2*kp)*L) + 2*Power(rm,2)*Power(rp,2)*Cos((km + 2*kp)*L) -        Power(rm,2)*Power(rp,4)*Cos((km + 2*kp)*L)))/((1 + 6*Power(rm,2) + Power(rm,4) - Power(-1 + Power(rm,2),2)*Cos(2*km*L))*(1 + 6*Power(rp,2) + Power(rp,4) - Power(-1 + Power(rp,2),2)*Cos(2*kp*L)))
		Ey = yr * math.cos(k0*z) - yi * math.sin(k0*z)
		z = zz
	else
		ar = -(((1 + rm)*(-Power(1 + rm,2) + Power(-1 + rm,2)*Cos(2*km*L)))/(Sqrt(2)*(1 + 6*Power(rm,2) + Power(rm,4) - Power(-1 + Power(rm,2),2)*Cos(2*km*L))))
		br = -(((1 + rp)*(-Power(1 + rp,2) + Power(-1 + rp,2)*Cos(2*kp*L)))/(Sqrt(2)*(1 + 6*Power(rp,2) + Power(rp,4) - Power(-1 + Power(rp,2),2)*Cos(2*kp*L))))
		cr = ((-1 + rm)*(-Power(-1 + rm,2) + Power(1 + rm,2)*Cos(2*km*L)))/(Sqrt(2)*(1 + 6*Power(rm,2) + Power(rm,4) - Power(-1 + Power(rm,2),2)*Cos(2*km*L)))
		dr = ((-1 + rp)*(-Power(-1 + rp,2) + Power(1 + rp,2)*Cos(2*kp*L)))/(Sqrt(2)*(1 + 6*Power(rp,2) + Power(rp,4) - Power(-1 + Power(rp,2),2)*Cos(2*kp*L)))
		ai = (Power(-1 + rm,2)*(1 + rm)*Sin(2*km*L))/(Sqrt(2)*(1 + 6*Power(rm,2) + Power(rm,4) - Power(-1 + Power(rm,2),2)*Cos(2*km*L)))
		bi = (Power(-1 + rp,2)*(1 + rp)*Sin(2*kp*L))/(Sqrt(2)*(1 + 6*Power(rp,2) + Power(rp,4) - Power(-1 + Power(rp,2),2)*Cos(2*kp*L)))
		ci = ((-1 + rm)*Power(1 + rm,2)*Sin(2*km*L))/(Sqrt(2)*(1 + 6*Power(rm,2) + Power(rm,4) - Power(-1 + Power(rm,2),2)*Cos(2*km*L)))
		di = ((-1 + rp)*Power(1 + rp,2)*Sin(2*kp*L))/(Sqrt(2)*(1 + 6*Power(rp,2) + Power(rp,4) - Power(-1 + Power(rp,2),2)*Cos(2*kp*L)))
		Ex = ar*math.cos(km*z)-ai*math.sin(km*z) + br*math.cos(kp*z)-bi*math.sin(kp*z) + cr*math.cos(-km*z)-ci*math.sin(-km*z) + dr*math.cos(-kp*z)-di*math.sin(-kp*z)
		Ey = ai*math.cos(km*z)+ar*math.sin(km*z) -(bi*math.cos(kp*z)+br*math.sin(kp*z))+ ci*math.cos(-km*z)+cr*math.sin(-km*z) -(di*math.cos(-kp*z)+dr*math.sin(-kp*z))
		Ex = Ex /math.sqrt(2)
		Ey = Ey /-math.sqrt(2)
	end
	
	print(z, Exr, Ex, Eyr, Ey)
end
